function [rebuiltModel] = rebuildModel(model,database,biomassReaction)
%
% Rebuilds a genome-scale reconstruction with Virtual Metabolic Human (VMH) 
% metabolic and reaction nomenclature while ensuring quality control through
% rBioNet.
%
% USAGE
% [rebuiltModel] = rebuildModel(model,database)
%
% INPUT
%    model            COBRA model structure
%    database         Structure containing rBioNet reaction and metabolite
%                     database
% OPTIONAL INPUT
%    biomassReaction  Biomass reaction abbreviation (if needs to be 
%                     specified, otherwise, will be inferred automatically)
% 
% OUTPUT
%    rebuiltModel  Quality-controlled COBRA model structure
%
% .. Authors
%       - Stefania Magnusdottir, 2016
%       - Almut Heinken, 12/2018: adapted to function.

% to account for older versions of AGORA
toReplace={'EX_4hpro(e)','EX_4hpro_LT(e)';'EX_indprp(e)','EX_ind3ppa(e)';'INDPRPt2r','IND3PPAt2r';'EX_adpcbl(e','EX_adocbl(e';'H202D','H2O2D'};

for i=1:size(toReplace,1)
    model.rxns=strrep(model.rxns,toReplace{i,1},toReplace{i,2});
end

% remove reactions not in database
notInDB=setdiff(model.rxns,database.reactions(:,1));
notInDB(find(strncmp(notInDB,'bio',3)),:)=[];
model=removeRxns(model,notInDB);

model=convertOldStyleModel(model);

% fix issues of incorrect field lengths
rxnFields={
    'grRules'
    'subSystems'
    'comments'
    'citations'
    'rxnECNumbers'
    'rxnKEGGID'
    'rxnConfidenceScores'
    };

for i=1:length(rxnFields)
if isfield(model,rxnFields{i}) && length(model.(rxnFields{i}))<length(model.rxns)
    if iscell(model.(rxnFields{i}))
    model.(rxnFields{i}){length(model.rxns),1}='';
    elseif isnumeric(model.(rxnFields{i}))
        model.(rxnFields{i})(length(model.rxns),1)=0;
    end
end
end
%% sort reactions
rbio=struct;
% get as much data as possible from the rBioNet database to avoid errors carrying over
for i=1:length(model.rxns)
    if ~strncmp('bio',model.rxns{i,1},3)
        % find reaction index
    rInd=find(ismember(database.reactions(:, 1), model.rxns{i,1}));
        model.rxns{i,1}=database.reactions{rInd, 1};
        model.grRules{i,1}=model.grRules{i};
        model.rxnNames{i,1}=database.reactions{rInd, 2};
        model.subSystems{i,1}=database.reactions{rInd, 11};
        if strcmp(database.reactions{rInd, 4},'1') || database.reactions{rInd, 4}==1
            model.lb(i,1)=-1000;
            model.ub(i,1)=1000;
        elseif strcmp(database.reactions{rInd, 4},'0') || database.reactions{rInd, 4}==0
            model.lb(i,1)=0;
            model.ub(i,1)=1000;
        end
        model.formulas{i,1}=database.reactions{rInd, 3};
    else
        model.rxns{i,1}=model.rxns{i};
        model.grRules{i,1}=model.grRules{i};
        model.rxnNames{i,1}=model.rxnNames{i};
        model.subSystems{i,1}='Biomass';
        model.lb(i,1)=model.lb(i);
        model.ub(i,1)=model.ub(i);
        model.formulas{i,1}=printRxnFormula(model,model.rxns{i,1});
    end
end

[uniqueRxns,oldInd]=unique(model.rxns);
rbio.data=cell(size(uniqueRxns,1),14);
rbio.data(:,1)=num2cell(ones(size(rbio.data,1),1));
rbio.data(:,2)=uniqueRxns;
rbio.data(:,3)=model.rxnNames(oldInd);
rbio.data(:,4)=model.formulas(oldInd);
rbio.data(:,6)=model.grRules(oldInd);
rbio.data(:,7)=num2cell(model.lb(oldInd));
rbio.data(:,8)=num2cell(model.ub(oldInd));
rbio.data(:,10)=model.subSystems(oldInd);
if isfield(model,'citations')
rbio.data(:,11)=model.citations(oldInd);
end
if isfield(model,'comments')
rbio.data(:,12)=model.comments(oldInd);
end
if isfield(model,'rxnECNumbers')
rbio.data(:,13)=model.rxnECNumbers(oldInd);
end
if isfield(model,'rxnKEGGID')
rbio.data(:,14)=model.rxnKEGGID(oldInd);
end
rbio.description=cell(7,1);

% build model with rBioNet
bInd = find(strncmp('bio',rbio.data(:,2),3));
if isempty(bInd)
    error('The model does not have a biomass objective function!')
end

for i=1:length(bInd)
    bAbb{i} = rbio.data{bInd(i),2};
    bForm{i} = rbio.data{bInd(i),4};
end
rbio.data(bInd,:) = [];
model = data2model(rbio.data,rbio.description,database);
for i=1:length(bInd)
    model = addReaction(model,bAbb{i},'reactionFormula',bForm{i}{1});%add translated biomass reaction
    model.comments{end+1,1} = '';
    model.citations{end+1,1} = '';
    model.rxnECNumbers{end+1,1} = '';
    model.rxnKEGGID{end+1,1} = '';
    model.rxnConfidenceScores{end+1,1} = '';
end

% fix incorrect format of PubChemID, metChEBIID, and metKEGGID
if isfield(model,'metPubChemID')
for i=1:length(model.metPubChemID)
model.metPubChemID{i,1}=char(string(model.metPubChemID{i,1}));
model.metChEBIID{i,1}=char(string(model.metChEBIID{i,1}));
model.metKEGGID{i,1}=char(string(model.metKEGGID{i,1}));
end
model.metPubChemID=cellstr(model.metPubChemID);
model.metChEBIID=cellstr(model.metChEBIID);
model.metKEGGID=cellstr(model.metKEGGID);
% fill in descriptions
model.description.author = 'Created by DEMETER, Molecular Systems Physiology group';
model.description.date=date;
end

% adapt model.subSystems to current standard
for i=1:length(model.subSystems)
    model.subSystems{i,1}=char(model.subSystems{i,1});
end

% set biomass reaction as objective function
if nargin >2
    model=changeObjective(model,biomassReaction);
else
    model=changeObjective(model,bAbb{1});
end

rebuiltModel=convertOldStyleModel(model);
end
    